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Abstract 


We present a practical cell-centred volume-of-fluid method developed within a pure Eulerian setting for the simula- 
tion of compressible solid-fluid problems. The method builds on a previously published diffuse-interface Godunov- 
type scheme through the addition of a specialised mixed-cell update that is capable of maintaining sharp interfaces 
indefinitely. The mixed-cell update is local and may be viewed as an interface-sharpening extension to the underly- 
ing diffuse-interface scheme along the lines of other techniques such as Tangent of Hyperbola INterface Capturing 
(THINC), and hence the method can be straightforwardly extended to include other coupled physics. We validate the 
method on a range of challenging test problems including a collapsing metal shell, cylinder impacts and the three- 
dimensional simulation of a buried explosive charge. Finally we demonstrate the robustness of the method, and its 
use in a multi-physics context, by modelling the BRL 105mm unconfined shaped charge with reactive high-explosive 
burn and rate-sensitive plasticity. 
© British Crown Owned Copyright 2023/AWE 
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1. Introduction 


Accurate and robust computational methods for the treatment of dynamic problems involving multiple solid and fluid 
materials are vital across many scientific and engineering disciplines. Applications include auto-mobile crashworthi- 
ness, sheet metal forming, inertial confinement fusion, asteroid impacts and supernova core-collapse to name but a 
few. Such problems frequently involve large deformations which cannot be treated by purely Lagrangian methods due 
to mesh tangling, and so evolving the interfaces between materials explicitly becomes a necessary part of any viable 
numerical approach. 

Many approaches have been developed to tackle this problem. Arbitrary Lagrangian-Eulerian (ALE) methods are 
designed to avoid mesh tangling by decoupling the motion of the computational mesh from the motion of the material. 
ALE codes are often implemented using a Lagrange-plus-remap approach where the mesh motion is formulated as 
an operator-split advection step appended to a pure Lagrangian update, but direct ALE is also possible [10, 4]. Fixed 
grid Eulerian methods on the other hand have many attractive characteristics, chief amongst them being the complete 
avoidance of all mesh-tangling related issues which can hamper robustness. Historically, many multi-material Eulerian 
methods have been developed based on solving hypoelastic systems on staggered grids, with explicit artificial viscosity 
terms to resolve shocks (see [10] for a review of such methods). In more recent years numerous Godunov-type 
methods, originally developed for compressible fluid flow, have been extended to solid dynamics using hyperelastic 
models [42, 36, 37, 31, 51, 8, 7, 18, 32, 5]. Godunov methods use the solution of Riemann problems to define 
numerical fluxes and introduce the required artificial viscosity implicitly. A distinct advantage of these pure Eulerian 
cell-centred methods is their amenability to implementation within adaptive mesh refinement (AMR) frameworks, 
which is often essential to render complex multi-physics problems computationally tractable. The foremost challenge 
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of integrating these methods into numerical schemes suitable for simulating multi-material problems is the treatment 
of material interfaces, which form internal boundaries. 

Methods for treating material interfaces can be broadly divided into two categories: interface-tracking and 
interface-capturing (analogously to methods for handling shocks). Interface-tracking methods use an indicator func- 
tion that denotes the location of the interface, such as the zero contour of a level-set [40]. Volume-of-fluid methods 
are another notable example, where the fraction of each material present in each cell is tracked and used to reconstruct 
an approximation of the interface [19, 41]. Cut cells may be discretised directly as polyhedra within a finite-volume 
or finite-element framework [37, 9], but this is complex and requires special treatments in order to avoid impractically 
small allowable timesteps. Alternatively, some approximate treatment may be used to evolve the solution around the 
interface, based on the true jump conditions. Examples include ghost-fluid methods, pioneered by Fedkiw et al. [17], 
where the state in each material is extrapolated across the interface (either by solution of the Riemann problem or 
some other method) allowing the use of an unmodified single-material update, which, by virtue of the manipulated 
states near the interface, captures the required interface boundary conditions [18, 6, 45, 32, 35]. A variety of closure 
models have been developed for performing Lagrangian updates in cut cells, typically based on asserting some kind 
of equilibrium between materials [46, 3, 24]. 

Interface-capturing methods differ in that the location of the interface is implicit in the solution. The update 
applied is the same irrespective of the number of materials present in each cell. The downside of this approach is 
that the interface tends to become smeared over time unless some sharpening procedure is used, hence these methods 
may also be referred to as diffuse-interface methods in contrast to the sharp (physical) interface usually offered by 
interface-tracking. A recent review can be found in [33]. Many diffuse-interface methods have been developed for 
multi-fluid systems [1, 44], but relatively few for solid-fluid interactions [16, 15, 5], principally because obtaining a 
consistent representation of the thermodynamic state in the mixture region is challenging due to the great variation in 
material properties. 

This paper is devoted to the development of a Godunov-type sharp-interface method that may be viewed as a hy- 
brid of interface-capturing and interface-tracking methods. Our method is based on the diffuse-interface scheme 
presented by Barton [5], which is itself a development of the five-equation multi-fluid model presented by Al- 
laire et al. [1] to include material strength. We design a specialised update procedure inspired by the multi-fluid 
methods of Miller and Puckett [38], and Cutforth et al. [12] for cells that contain more than one material, which 
preserves the sharp character of interfaces by means of volume-of-fluid advection. The solution within these cells is 
derived from the Riemann problem posed in terms of the underlying diffuse-interface model, and a pressure relaxation 
step is included to force the solution to converge to this model. Our mixed-cell update may be viewed as an alter- 
native to other interface sharpening methods such as Tangent of Hyperbola INterface Capturing (THINC) [60, 49] 
(previously applied to this problem by Barton [5]) and other anti-diffusion schemes [48, 50, 25]. The method is prac- 
tical (not requiring any complicated geometric considerations), robust (capable of treating multiple solids and fluids 
undergoing extreme deformation), and extensible (amenable to the addition of coupled physics). 

The remainder of the paper is organised as follows. In Section 2 we lay out the theoretical underpinnings of 
the method, including the evolution equations and the thermodynamic and kinematic framework. In Section 3 we 
describe the numerical method in detail. In Section 4 we validate the method through application to a strenuous series 
of problems culminating in the simulation of an explosive charge buried in clay. Finally in Section 5 we show the 
potential of our method by using it to model the BRL 105mm unconfined shaped charge, an extremely challenging 
multi-physics problem featuring reactive high-explosive burn and rate-sensitive elasto-plastic solid dynamics. We 
conclude the work in Section 6. 


2. Governing theory 


2.1. Evolution equations—derivation 


The sharp-interface method is presented as an extension to the diffuse-interface method of Barton [5]. Under this 
model, materials are allowed to mix at their interfaces. Each material / is described by its mass density p', volume- 
fraction ¢! and specific internal energy &’. The model assumes mechanical equilibrium: materials in a mixture share 
a common velocity u and mixture stress tensor o. The basic conservation equations for mass, momentum and energy 
are then as follows 
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a TV (pu) = 0 ” 
OWE) (peu) = V- (ow) 3) 


where p, pDE = p& + p(u-u)/2 and & are the mixture density, total energy density and specific internal energy 
respectively, which along with the volume-fraction are subject to the thermodynamically consistent mixture rules 


1 := Sie! (4) 
1 
p = dp! (5) 
X 
pe = Sole! (6) 
1 


In contrast to the model from Barton [5] which uses simple advection equations for volume-fraction, here we use 
the volume-fraction evolution equation from [38, 31, 12] which is derived directly from the continuity equation (1) 
and therefore takes account of the relative compressibility of each material due to the divergence of the velocity field. 
First, expanding product derivatives gives 


og! I I 

se MGW ep (7) 
If we assume that any compression that takes place is isentropic, and that isotropic stress is maintained during the 
advection process, then the pressure change in each material will be equal to the pressure change of the mixture. Then 
we can use the definition of the isentropic bulk modulus Ks = p 0p/0p|, to replace references to the partial densities 


in (7) with fractions of the mixture density. 


K 
jeg, (8) 
ve) 


A 41 ly a 
Be Gay |e +a. v9] (9) 


Finally summing (1) over materials to give the continuity equation of the mixture and substituting 0p/0t + u- Vo = 
—pV -u into (9), we obtain the volume-fraction evolution equation! 


og! 1, _ &Ks 
ar + ¥: (eu) = xi V-u (10) 


'We note that while the choice of right-hand side weighting in the volume-fraction equation can be motivated physically as described, phe- 
nomenological choices can also work well. We have tested substituting the Griineisen parameter for the isentropic bulk modulus, to make the 
volume-fraction evolution appear consistent with the internal energy evolution described later, and found no significant difference in the perfor- 
mance of the method. 
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The pressure for each material is assumed to be defined by an equation of state of the Mie-Griineisen form 


p'(o!, &',dev(H,)) = Pree(0', dev(H)) + p'T'(p')(! — Exee(0!, dev (H2))) (11) 


where I” is the Griineisen parameter and dev(H,) is the deviatoric part? of the 3 x 3 Hencky elastic strain tensor. The 
reference energy is assumed to admit an additive decomposition, comprising a contribution due to cold compression 
or dilation, &, and a contribution due to isochoric shear strain, &, 


Erer(p',dev(He)) = &(o') + &(p!, dev(H.)) (12) 
lol l 01 OS ret 
Pret(P -dev(H,)) = pay (13) 
cp 
The shear energy is given by 
Il l G'(p') 2 i 
&5(p', dev(H,)) = ose (dev(H.)) (14) 


where G! is the shear modulus and .J*(dev(H,)) = tr(dev(H,)dev(H,)’) is the second invariant of shear strain. In 
fluids the shear modulus is taken to be zero, leading to zero shear energy contribution. 
Following hyperelastic theory the Cauchy stress tensor is a product of thermodynamic compatibility and is found 
through taking derivatives of internal energy with respect to density and deviatoric strain (see [5] for a full derivation) 
o!(p', &', dev(H!)) = —p!(p!, &', dev(H!))I + 2G'(p')dev(H) (15) 
Using the assumption that all components in a mixture are in mechanical equilibrium and therefore have the same 
pressure, along with the assumption that materials have a common deviatoric strain, the mixture stress becomes 


a(p, &, dev(H.)) = —p(p, &, dev(H_.))I + 2G(p)dev(H. ) (16) 


where the mixture pressure p and mixture shear modulus G are given by 


p& —S (¢'p'El, — o'pl,/T') 


= (17) 
we /T! 
ig! TV’ 
Gs 2G /P (18) 
di e/T 

The mixture Griineisen parameter I is defined as 

= 1 

= (19) 


¢/T! 


In order to track elastic deformations in solid media, it is first noted that the Hencky deviatoric strain can be 
defined in terms of the unimodular elastic left stretch tensor V, according to 


>The matrix deviator dev(M) is defined as M — 4tr(M)I, where tr(M) denotes the matrix trace. 
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dev(H,) = In(V,) (20) 


An additional nine equations are then included to evolve V. 


AW = 2 _ 
ONE 1 a) NaN s(V-u)V,—® (21) 


The source term ® represents relaxation due to inelastic deformation, defined later. 

It is sometimes necessary to evolve additional history variables to support the closure models, For instance equiv- 
alent plastic strain or the reactant mass fraction in a condensed-phase explosive. These additional history parameters 
can be collectively denoted by the vectors ! for each material, and evolve according to 


ue 


+V. (p'¢'x! @ u) = pigla! (22) 


For the interested reader, further examples of how multi-physics can be introduced into the model in this way include 
the treatment of reactive mixtures in [55], continuum damage mechanics and fracture in [56], and phase transitions 
in [59]. 

Finally, around material interfaces we require an additional set of equations to evolve the partial internal energies 
of each material. In [38] an internal energy evolution equation was proposed which added a pdv work term to the 
right-hand-side of the conservative advection equation, inspired by the right-hand-side of the volume-fraction equation 
(10). Here we propose an alternative formulation which is equally self-consistent, but which uses the mixture rules to 
determine the additional source terms, including those taking account of material strength. By expanding the energy 
equation (3) into internal and kinetic contributions we obtain 


(PE) 
Ot 


+V-(p&u) —o0 : Vu=0 (23) 


Using the definition of stress (16) and substituting the mixture rules for energy, pressure and the shear modulus gives 


gT 
u) — ae! + 2G'dev(H.)) : Vu} =0 (24) 


from which we assume the internal energy equations for individual materials 


a al nL el 
Ce 


- —p'l + 2G'dev(H.)) : Vu (25) 
0 


gT 
V- (¢'p'é"u) = oa 
2.2. Evolution equations—summary 


The complete system of evolution equations for which we seek solutions comprises (1)-(3), (10), (21), and (22). In 
vector form, this system can be written: 


oU OF; ( 
1p tier, a = LSU ) + Syptie(U) (26) 


Ot 


where the vector of conserved variables U and the fluxes in the x, direction F,(U) are defined (in Cartesian coordi- 
nates) as 
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pg! p' plu 


¢! ug 
_ pu a PUujug — ix 
CS pE |? F,(U) = pEuy — >); Citi (27) 
Vei i Veji je a Veo gjUi 
p v) Tin Pp wv) Ty Uk 
The source terms are separated reflecting the splitting strategy used in the numerical implementation 
0 0 
¢'Ks/Ks - Oug/ OX 0 
0 0 
Sx(U) = 0 > Ssi(U) =| 4 (28) 
2V.ijOun/OXE = UiOV ex /OXk Sia 
0 POTn 


To assist in determining solutions in mixed cells, the system is complemented by the evolution equations for internal 
energies (25). 


2.3. Quasi-linear primitive formulation and wave speeds 


Our choice of numerical method requires solution of the quasi-linear primitive form of the evolution equations, which 
neglecting the inelastic and history-rate source terms can be written 


OW OW 
+ > A.(W)— =0 29 
a LNW), (29) 


where W is the vector of primitive variables 


W = (0'6'».6'.ui, D, Ves Tn) (30) 


The matrix in the quasi-linear formulation, A, is given by (in the x-direction) 


vu 0 pe 0 0 00000000000 
Ou 0 0 0 0 0 0 0 0 0 0 0 0 00 
BY Cy u 0 0 Ip Ay, Ay Ady Aly Ayy Ajy Aj Aj3 A33 0 
By C; 0 ‘eg 0 0 At, A3, Az, Aj, AS) AZ, Aj, 43, A, 0 
BC} 0 0 " 0 A}, A}, A3, Aly Ady A3y A}; 43, 43, 0 
0 0 pa 0 0 u 0 0 0 0 0 0 0 0 00 
00-21 0 0 0 « 0000000 00 
00 ¥en -Ver 0 0 0 u 0 00000 00 
A.(W) =] 00 '¥%5 0 —Kn 0 00 «0000000 (31) 
00-23%» 0 0 0000 u 0000 00 
00 tv» -Ver 0 0 000 0 u 00 0 00 
00 1% 0 -Var 0 00000 u 00 00 
00-2%53 0 0 0000000 u 0 00 
00 4¥e3 -Ver 0 0 0000000 u 00 
00 1¥33 0 -Va3 0 00000000 u0 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 04 
where 
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a Odev(a);1 


C= = (32) 
BO OB GV 
1 ddev(a); 
pi = 1 Movie (33) 
p Op'¢') 
1 od ; 
gen (34) 
p og! 
(35) 
The derivative of the tensor logarithm is given by Jog [21], and so 
—-1 
Ai, = —2b°6ij5uVein (36) 


The bulk sound speed a and the shear speed b are defined next. 

The wave speeds are required by the primitive variable formulation, the Riemann solver, and in order to estimate 
the allowable time-step. In [5] it was shown how the non-linear characteristic wave speeds derived from the charac- 
teristic polynomial of A, are a function of the eigenvalues of the 3 x 3 acoustic tensor. In lieu of a costly evaluation 
of this tensor the fastest longitudinal wave speed is estimated as 


4 
c"(p!, p', dev(H,)) = a*"(p!, p', dev(HL)) + ae (p') (37) 


The bulk sound speed is computed by differentiating (11) with respect to density at constant entropy, and using the 
relation d& = TdS — pdv 


op! ,(p', dev(H,)) : (I'(p') + 1)p' — pl (e', dev(H.)) | 


a’(p!, p',dev(H,)) = (38) 
ap! p! 
! = Prego! dev(H)) OF ey | 1 OSee( dev HH) 
Mp!) foe RY ap! 
The shear speed is given by 
ppl) — EL) (39) 
~ 1 
p 
In cells containing more than one material, the mixture wave speed is calculated according to 
Lol 2p! 
2 — whee fe (40) 


yt 


2.4. Closure models 

It is noted that the Mie-Griineisen framework embodies many specific equations of state, depending on the functional 
form assumed for the reference energy/pressure, Griineisen parameter and shear modulus. For example (superscripts 
l have been omitted for clarity): 


e The ideal gas law is given by él. = Pret = 0, =y—1landG =0. 


e The stiffened gas equation of state describing denser fluids is given by Step = Ec, Pref = —YPo,T = y — 1 and 
G=0. 
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e The Jones-Wilkins-Lee (JWL) equation of state [27, 29, 34] describing explosive products (and sometimes 
reactants) is given by [ = 9, G = 0 and 


A B 
Gref = exp ( ni) + exp ( rf) 
Ripo p Rp p 


e The equation of state from [13] describing elastic solids is given by T = I and 


a 2 
K 
Eref = = ((4) 1) +6; 
2pow Po 


> Oref 
Pref = P Tae 
Ye) 


B+1 
o-0( 5) 
Po 


Plastic relaxations are incorporated via the source term ® following the method of convex potentials. The Von 
Mises flow rule is used, leading to the following functional form 


3 devia) = 
V. 41 
2 [devo re 


D=y 


where ||M|| = (M) is the Frobenius norm and y > 0 is the plastic rate. The latter is a closure model which must be 
provided for each material, and the following mixture rule is used 


pe py gy /T! 
a e/T! 


In this work we consider both ideal plasticity, and the rate-dependent model from Johnson and Cook [23]. For the 
former the plastic rate is given by 


3 
x! = XoH ( ne IIdev(o)|| — n) (43) 


where H(s) is the Heaviside function, oy is the constant yield stress, and x — oo. For the rate-dependent model, we 
have 


(42) 


1 a 3 ||dev(o)|| _ 
X = Xoexp E ( ey \) (44) 


where xX is the reference rate, c3 controls the rate dependency, c, is the yield stress, c is the strain hardening factor, 
n is the strain hardening exponent, Tyner is the melting temperature of the material, To is a reference temperature, and 
m is the thermal softening exponent. This model is a function of the equivalent plastic strain Eb so an additional 
evolution equation is required to track this quantity 
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og'p'e’, 
ot 
which is added as a history variable for the corresponding material. 
Condensed-phase explosives are treated using the reactive burn model from [55], where the reactants and products 


are represented as a physical mixture. The reactant mass fraction, denoted A, is added as a history variable for the 
corresponding material with the evolution equation 


+V-(d'p'e,u) = ¢'p'y! (46) 


a¢'pla! 
Ot 


+V-(¢plalu) = gpla! (47) 


Various forms for the reaction rate have been proposed, notably the ignition and growth model [30], viz. 


jen (2 ae a) H(Fg —f) + (48) 
Po 

G,A° f'p (Fo, — f) + (49) 

G2A° fe pA f — Fa) (50) 

f=1Su (51) 


where H(s) denotes the Heaviside function. 


3. Numerical method 


We leverage the AMReX block-structured adaptive mesh refinement framework in our implementation [14]. AMReX 
has native support for advanced architectures such as Graphics Processing Units (GPUs) which have the potential to 
dramatically accelerate calculations. All of the 2D and 3D simulations presented later in this work are performed on 
NVIDIA A100 GPUs. 

Solutions are found using structured computational grids consisting of hexahedral cells with cell centres denoted 
by the indices i, j,k. Each cell Cj; has the dimensions Ax; = Xj41/2 — Xi-1/2, AYj = Yj41/2 — Yj-1/2, A& = 
Zk41/2 — Zk-1/2» Where i + 5, ifac ‘, k+ 5 denote cell boundary quantities. Each cell therefore forms the control 
volume V;;, = Ax;Ay;Az,. The method is implemented in a dimensionally split fashion, which ensures coupling at 
cell corners is accounted for, and also greatly simplifies the volume-of-fluid advection. A timestep thus consists of x, 
y and z updates, followed by the split source term updates. Within each cell the solution is advanced according to 


n(1 Ere At ii 
u, _ Uin - AG, [Fiitt/2 — Fii-12] + AS (UF) 
n(2) _ yyn(t) _ At n,(1) 
Ui ~ Uri ~ Ay; [Fo j+1/2 ~ F» j-1/2| ci AtS2(U; 5, ) 
n(3) _ zqn(2) _ At »,(2) 
Ui = Ui, ~ Az, [Fsu+1/2 a F34-1/2| i ArS3(U; i, ) 
Ung! = Uj) + ASsou(Ujn”) 


where F; ;+1/2 are the cell-wall numerical fluxes to be defined later. The volume-of-fluid method is incoporated into 
the dimensionally split part of the update, described in detail below. 
The basic procedure may be outlined as follows: 


1. For each coordinate direction: 
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(a) Perform MUSCL reconstruction of the primitive variables and extrapolate to the half-time-step according 
to the MUSCL-Hancock method. 

(b) Use the HLLD Riemann solver [32, 5] to compute interfacial states and numerical fluxes at all faces. 

(c) Reconstruct material interfaces using the initial volume-fraction field and use the interfacial speeds from 
step 1b to compute the volume of each material transported across each face. 

(d) Local to material interfaces, apply the mixed-cell update formulae to evolve the state. Then relax each 
mixed-cell to pressure equilibrium and sum the partial internal energies to retrieve the total energy. 

(e) Update state within single-material cells (i.e. all cells not updated in step 1d) using the HLLD fluxes from 
step 1b. 


2. Calculate the operator-split source terms. Reinstate the symmetry property for the unimodular elastic left stretch 
tensor. 


For clarity we describe the update in the x-direction, modifications for the y- and z-directions are trivial. 


3.1. Reconstruction 


In order to obtain second-order accuracy in smooth single-material regions, we use the MUSCL-Hancock method [53], 
wherein the input data for the Riemann problems is computed by performing a limited MUSCL reconstruction of the 
initial cell average states, and then extrapolating the reconstructed states forward in time to the half-step ¢ + Ar/2. 
These states are then used to compute the numerical fluxes. The method can be applied to either the conserved or 
primitive variables, but testing has shown that reconstruction of the primitive variables yields better results, with fewer 
oscillations in the vicinity of material interfaces. This is in line with the findings of Johnsen and Colonius [22]. 

The reconstruction is carried out as follows. Given the piecewise constant vector of primitive variables W in each 
cell, the MUSCL method constructs a piecewise linear function from which the left and right boundary values can be 
written 


1 


Wi-i2r = Wi- 5 P(Ar) (Wi41 — Wi), (52) 
1 
Wi+1/2,1 — W; + 5 P(Az) (W; = W,-1) > (53) 
with 
Wi+1 = W; Wi _ Wi-1 
A, = =———_. = 54 
E W; — Wi-i + € m Wi4i —Wit+e (54) 


where € is a vanishingly small number used to avoid division by zero. The function ®(A) represents the slope limiter, 
which prevents spurious oscillations and ensures that the Total Variation Diminishing (TVD) property is maintained. 
In the examples presented here the limiter of Van-Leer [52] is used 


_ A+ |A| 


The extrapolation in time is carried out using 


n+1/2 n At n n n 

Wier = Wier 7 Aga (Wi) [w: = We pal (56) 
n+1/2 _ n At n n n 

Wien = Wiss. — Ah (Wi) hare = w;| (57) 


We note that in [5] material interfaces were sharpened following MUSCL reconstruction using the algebraic 
Tangent of Hyperbola INterface Capturing (THINC) technique from Xiao et al. [60], but that this is unnecessary in our 
scheme as the interfaces are always confined to a single cell by virtue of the specialised mixed-cell update. However 
it may still be advantageous to apply the THINC technique to fields other than the material volume-fractions. For 
example, Wallis et al. [56] use THINC to sharpen the material damage field used for Continuum Damage Modelling 
(CDM). 
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Figure 1: The HLLD solver assumes the solution to the Riemann problem contains six constant states interposed by five discontinuous waves: two 
fast longitudinal waves, two slow shear waves, and one central contact discontinuity. 


3.2. Solution of Riemann problems 


The numerical fluxes and interfacial states are obtained using the HLLD approximate Riemann solver, originally 
developed for magnetohydrodynamics but repurposed for solid dynamics in [32, 5]. The solution consists of six 


1 
constant states interposed by five waves (see Figure 1). Given the reconstructed left and right states U, (wr >’) and 


1 
Ur(Wi, >) as input, the fluxes are 


F(U,) if: 525.0 
F(U,) + S, (Uz — U;) oS. 20<5* 
F(U;) + Si (Us — Uz) + S3}(U7*—Us) if Si <0<S* 
Fuiip = E et Sere ‘ i : (58) 
F(Ug) + Sp(UR — Ur) + S,(UR* — UR) if S*<0<S5 
F(Up) + S(Ug — Ur) if Si<0<s} 
F (Up) if Sh 0 


where the left longitudinal, left shear, contact, right shear, and right longitudinal wave speeds are given respectively 
by 


Ss a min(uz, — Clr, UR — CR) (59) 
$i =s* bt (60) 
uz(S! —uz) — prur(S4 — ur) +or11—-o 
s* _ PL L( E Zl PR R( R eo Ll R11 (61) 
pPL(S7, — UL) — pr(S'p — UR) 

Si =S* +3 ae 

Se = max(uz + Cr,UR + Cr) (63) 
Note that in fluids, b = 0 and therefore Sj; = S* = S; so the solver reduces to the standard HLLC scheme. 
The intermediate states U7, U7*, U;* and Uj, are obtained through manipulation of the jump conditions 

Si (Ux — U%) = F(Ux) — F(U%) (64) 

Si(Ux — UX") = F(Ux) — F(UZ’) (65) 


It is assumed that normal stress and normal velocity are constant across shear waves, but may jump across longitudinal 
waves. Similarly, the transverse stresses and velocities may jump across shear waves but are assumed constant across 
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longitudinal waves. Solving (64) under these assumptions gives the states between the two pairs of longitudinal and 
shear waves 


PRWK 
pEK + (S* = UK) (pKS* = oxi/(Se = ux)) 
VeK/XK 
: VeK21 
Uk =XxK Viera (66) 
Ve,K.12/XK 
Ve,K,22 
_ Ve,k,32 
Ve,K,13/XK 
Ve,K.23 


Ve.K33 
Laid 
PRO KT Km 


where 


Si. —uK 


= oh (67) 
si=s 


XK 


Between the two shear waves no-slip boundary conditions—constant transverse stresses and transverse velocities— 
are used to describe how the solution varies across the contact wave, which reflects the underpinning mechanical 
equilibrium of the model. The intermediate states then become 


0 
0 
0 
OK, — TK21 
ORs) — TK 31 
VE OKI —VKOK.1 + WOR 31 —WkOK31 


0 
1 Ve cuve* — vx) 
=U ee 68 
x SP Se Vex — WK) 
Vex io(Yx — vx) 
Veni(we — wx) 


0) 
—* 
Vex e — vx) 
Veniwe — wx) 
0 


We also require the intermediate stresses which are given by 
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Ole = TUK — PK(S — UK)(S* — ux) (69) 


Cnn =x for j = 2,3 (70) 

OK = TiLK (71) 

ae win) Se SE ep 8, 3 (72) 
where 


ax = pk(S* — 8%) 


The fluxes Fup are used directly to define the cell wall fluxes F;;+1/2 to update single-material cells following 
the standard MUSCL-Hancock method. In order to perform the mixed-cell update the fluxes are not needed, but we 
do require the states evaluated at the interface 


U, if si>o 
Uy a S057 
UF if St <0<5* 
UPS ea: (73) 
U, if S*<0<S, 
Us. if Si<0< 8) 
Ur if 5 <0 


It is important to note that where the interfacial stresses are required (e.g. in (84)), these must be taken directly from 
(69)-(72). Evaluating them from U™ using the equation of state is invalid and will cause the scheme to fail. 


3.3. Interface reconstruction 


Material interfaces are reconstructed using the method of Youngs [61, 41]. In each cell, and for each material, the 
gradient of the volume fraction field is approximated using finite differences and an outward-facing normal vector is 
calculated 


1 
Yee (74) 

IV¢"| 
An oriented plane associated with this normal is then positioned such that it cuts the cell producing the required 
volume fraction. This plane is taken as a linear approximation of the interface within the cell. In cells with two 
materials, the second material interface is simply the first with the sign of the normal vector reversed. In cells with 
three or more materials, producing an accurate reconstruction is more complicated and has been studied by many 
authors [11, 26, 47]. This problem is outside the scope of the current work, and we adopt a simple approach. The 
user of the code specifies a global ordering of materials, and interfaces are computed between adjacent pairs in this 
ordering using a cumulative sum of the individual volume fraction fields. For instance, with four materials, we would 
compute the interface of material one using ¢!, material two using ¢! + ¢? and material three using ¢' + ¢? + ¢°. The 

interface of material four would be the same as material three with the normal sign reversed. 

Once we have the required material interfaces defined, we require the signed transported volumes of each mate- 
rial across the left and right faces of the cell (in the sweep direction), AVI 12 and Av! +1/2 respectively. These are 
computed as the intersection of the material polyhedra underneath the interfaces and the cuboidal region traced out by 
the motion of the continuum normal to and through the cell face by a distance Ar u'", with sign equivalent to wi". In 
the x-direction this cuboid has the signed volume AyAzArui™, and the sum of the transported volumes 57, AV! must 
exactly equal this value at each face. A 2D illustration of this process is given in Figure 2. 
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Figure 2: Illustration of the calculation of transported material volumes at each face. The filled regions show the reconstructed interface between 
two materials. The hatched regions denote the advected volumes, with widths given in terms of the interfacial speeds. The transported material 
volumes consist of the intersection of these regions with each material. 


3.4. Mixed-cell update 


The specialised mixed-cell update must be applied in all cells containing more than one material, but additionally in 
cells that could become mixed following the update. In practice this means all pure cells that are adjacent to a mixed- 
cell in the direction of the sweep, or two adjacent pure cells containing different materials (in the case where the 
material interface coincides exactly with the adjoining face). Even in the most demanding multi-material calculations 
this will typically be less than 10% of the total cell count, so the computational burden of this process is small. 

Following interface reconstruction, the intermediate volume fractions (neglecting contributions from the right- 
hand side of (10)) are computed using volume-of-fluid advection 


AV; — AV; 
V 


# =!" (75) 


where V is the volume of the cell. In this section the subscripts L and R are used as short-hand for i — 5 and i + 5 
respectively to denote quantities located at the left and right faces of the cell under consideration. 

Following application of (75), the sum of the volume fractions ¢! in the cell may not equal one (whenever the 
velocity divergence is non-zero). The difference (1 — 5}, 6”) must be redistributed amongst the materials in the cell 
in accordance with their relative compressibility. To compute the right-hand side of (10) we require the bulk modulus 
for each material following advection, denoted K!., which we obtain through a volume-weighted averaging procedure 
as detailed in [38, 12] 


AV KS, + Gi VIKS a AVER 


gl 

7 76 
S,i EV; ( ) 
where the upwind bulk moduli at the left and right faces are given by 

a Ke if act 0 

Sr a m2 (77) 
, Ky F otherwise 

| K ‘ F if wv, 0 

Ky r= ; +3 (78) 
, K : i, otherwise 


The post-advection bulk modulus of the mixture is then 
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-1 
- | 
Kee BaP (79) 
(=f 


Note that the same approach is used to determine the post-advection Griineisen parameters I” and I used in the internal 
energy equations below. The final volume fractions are then given by 


‘ as EK 
gmt = Fa (1 gS = 
m S 


TM\ __ UR Uy, 
(1 28 )= AA (81) 


The mass, momentum and energy transported across the face of each mixed-cell are computed using upwinding. 
This is important to ensure that the procedure is self-consistent: when the cell is emptied of material /, the correspond- 
ing contributions to mass, momentum and energy must be set to zero. The relevant upwind quantities (denoted with 
a wide hat q) are calculated at the left and right faces of the cell as in (77) and (78) respectively. We also require the 
transported density of the mixture at each face, given by 


~l ayl 

oo= Du PKAV 

Ke pe 
dU AVE 

The discrete update equations for the partial densities, momenta, elastic stretches and history variables are as 


follows. The advective terms are differenced in line with (75). Note that the inelastic and history-rate source terms 
are omitted as these are evaluated separately. 


, forK=L,R (82) 


Avie ~ Avia 


occ pee = pin 7 (83) 
AV aprile — AV Pit Cir ~ Tit 
+1, n4+1 RPRULR LPLYLL i1,R il,L 
pat = ptt — a — 
a AVaVeiir -_ AVL Veit 2—n 41 
Ve = Ver V 3 ei j(l S38) oa (85) 
1 
At . int . int 
ae [Oe — wi Veie + (ute — 0)Ver | 
InI4I IAI AI 
Girt lint glntl — glnplngln AVRPRT nk ~ AVP mt (86) 


m V 


The momentum update (84) is importantly not the same as the one from [38], which used the interfacial velocity rather 
than the upwinded velocity in the advective terms. We found the original to be highly unstable when a dense material 
is leaving a cell containing a light material, as the dense material’s contribution to the momentum is not correctly set 
to zero. This resulted in large spikes in velocity which caused calculations to fail. The new update does not exhibit 
these problems. 

The discrete energy update based on (25) is given by 
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geipea = pple tr? Aadv = pe Avei + ee (87) 
avigdl — avid 


Aa v= 
; V 
7m oT 
Avo = (1 »? ) fl 
Ik re ume — un | éT 
Aghear = 2G” d dev(H.) 7; ge ar 


Cutforth et al. presented an energy update which reportedly improved robustness significantly in complex simula- 
tions [12]. Following their approach, the time-level of the pressure p'* in the discrete energy equation is taken to 
differ depending on whether the cell is in compression or expansion: 


in . int int 
pit = . if Sy ay (88) 


p'"*! — otherwise 


In order to ensure a consistent evaluation of the stress, the same logic is applied when chosing the time-level of the 
shear modulus G!* and the deviatoric strain dev(H,)*. 

In the case of compression the update may be readily evaluated. In the implicit case (expansion), we invoke the 
equation of state (11) to rearrange the right-hand side 


Ln olan gl, Lat pint glint Int 
pp" — Baa + Avoi(e gaa Beas Err —p ) + Aghear 


Intl goln +1 ref 89 
T+ Ayo lt! /plint! ( ) 


pirtly 


Then all of the variables on the right-hand side of (89) are known: ¢/"+! is given by (80), p!”*! is given by (83), 
dev(H.)"*! is given by (85), and T"*!, ets po and G!"*! are all assumed to depend only on p!”+! and 
dev(H,)"*". 

In the case of other equations of state where a closed-form rearrangement is not possible (for instance reactive 
mixtures along the lines of [55], or cavitating liquids along the lines of [59]), a root-finding procedure may be used to 
perform the implicit update. 


= 


3.5. Pressure relaxation 
After the update mixed-cells will not in general be in pressure equilibrium, which is a requirement of the base diffuse- 
interface method. A pressure relaxation step is therefore included to iteratively adjust volume-fractions and internal 
energies until equilibrium is attained. We first outline the pressure relaxation scheme proposed by Miller and Puckett 
in [38], and then describe changes we have made to improve robustness. 

The following equations are solved to adjust the component pressures in mixed-cells towards equilibrium in a 
controlled manner. 


Ky 

of 

Sj Ag! = 0 (91) 
1 


p= p'+Ap'=p'-—A¢' (90) 


These yield approximate expressions for the equilibrium pressure p and the required change in volume-fraction for 
each material Ag! 
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fi om, ¢'p'/ Ky (92) 
mY 
an 
Ag! = 57 (P'—P) (93) 
S 
Due to the linearisation implicit in (93), Miller and Puckett enforced a limit on the maximum relative volume change 
in a single iteration of the relaxation, viz. |Ad'/¢'| < 6. When the cell is under compression 6 = 0.1, and 6 = 0.05 in 
expansion. 
On each iteration, (92) is used to obtain an estimate of the equilibrium pressure, which is then used in (93) to 
calculate the amount by which each volume-fraction should change. If any of the values exceed the threshold, all of 


them are rescaled (this ensures that they continue to sum to zero) 


6 
Ad! — Ad! min ——— (94) 
1 |A¢'/¢"| 
Finally the volume-fractions are changed, and a pdv correction is made to each partial internal energy 
og + Ag (95) 
pie! — p'glé! — pad! (96) 


The process is repeated until each component pressure is within some prescribed tolerance. 

We found the above scheme to lack robustness, particularly within solids where tiny changes in volume-fraction 
result in relatively large changes in pressure. We observed oscillatory behaviour around the desired equilibrium 
point, where the calculated volume-fraction changes would repeatedly overshoot the target pressure, resulting in non- 
convergence. In order to address this we modify the threshold adaptively. Each time a sign-change is detected in Ad! 
we reduce the threshold value 6. The threshold is given by 


5 = doexp(—aNgign) (97) 


where 6o is the initial threshold value, a is a free parameter that controls the rate of fall-off, and Ngign is the number of 
iterations in which a sign-change in Ad! has been detected. We typically take 59 = 0.01 and a = 0.1. Unlike Miller 
and Puckett we do not distinguish between compression and expansion when applying the threshold. 


3.6. Plastic update 


The plastic source term ® is evaluated as detailed by Barton [5]. This update is independent of the volume-of-fluid 
method and is performed identically in pure and mixed cells, but is summarised here for completeness. First the 
following ordinary differential equation is solved using the backwards Euler method to evolve the strain invariant for 
each material that adheres to a plasticity law 


. 3 
ZS! 225 [3x gs” E (0, 7°] (98) 


and the mixture strain invariant is obtained following 


ee 
o> gt 
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(99) 


Material po (g cm~3) Ky (GPa) Go (GPa) To a p 


Aluminium 6061-T6 2.703 76.3 26.36 1.484 0.627 2.288 
Copper 8.930 136.45 39.38 2.0 1.0 3.0 
Clay 1.932 49.4 6.0 1.0 1.0 - 


Table 1: Material parameters for the equation of state from [13], as given in Section 2.4. Where £ is not specified a constant shear modulus is used: 
G =Go. 


where J! = 0 in fluids. It is noted that since the ODE only need be solved for materials where ¢! > 0 the volume- 
of-fluid method increases the efficiency of this update by virtue of there being fewer mixed cells. Finally the mixture 
strain invariant at the new time level is used to update the stretch tensor 


V, = exp (= dev(H,) (100) 


In the course of this last step the symmetry and unimodularity of the stretch tensor is reinstated, which may be lost 
during the hyperbolic update: 


V.<— VV. (101) 
4. Verification and validation 


4.1. Solid-solid Riemann problem 


The first test is a solid-solid Riemann problem taken from [8], involving an interaction between aluminium and copper 
resulting in a complex wave structure. The initial left and right states are as follows: 


2 1 0 0 
u,=[ 0 |kms~!, F.p=[ -001 0.95 0.02), & = Getz 
0.1 -0.015 0 09 
0 1 0 0 
ug = | -0.03 }kms~!, Fire= {0.015 0.95 0], &= Gere 
—0.01 -0.01 0 09 


where F, is the elastic deformation gradient tensor, from which the initial densities and elastic stretches are found via 


Po 


2 102 

p~ det|F.| (02) 

_ 4 /E.FT 

V,= 5 (103) 
\/det|F.| 


Both materials are assumed to be purely elastic and are governed by the equation of state from [13] with parameters 
given in Table 1. The domain is x € [0, 1] with the discontinuity initially located at x = 0.5, with aluminium to the 
left and copper to the right. The problem is run to time 0.5 zs with CFL number 0.1. 

Results are presented in Figure 3 for a mesh with 500 cells. The volume-of-fluid method (labelled VOF) is 
compared against the base diffuse-interface method with THINC interface sharpening [5] (labelled DI+THINC) and 
the exact solution derived using the method from [8]. Convergence is illustrated in Figure 4, which gives L,-norms 
for density and velocity relative to the exact solution. 

As expected, away from the material interface the results are equivalent to the diffuse-interface solution, however 
at the contact the ability of the volume-of-fluid scheme to maintain a sharp interface is demonstrated. Barton noted 
that THINC contributes to oscillations in the tangential stress and velocity fields at the interface—these are notably 
reduced with the new method. 
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Figure 3: Solid-solid Riemann problem. 
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Figure 3: Solid-solid Riemann problem (cont.) 
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Figure 4: Solid-solid Riemann problem, L;-norms versus exact solution. 
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4.2. Solid-fluid Riemann problem 

The second test is a solid-fluid Riemann problem, again taken from [8], featuring highly stressed moving aluminium 
in contact with quiescent air. This test is used to verify that the volume-of-fluid method can correctly treat an interface 
between two materials with radically different properties. The initial left and right states are as follows: 


2 1 0 0 
uz=| 0 |kms-, F.p= | —001 0.95 0.02], & = Getz 
0.1 -0.015 0 09 


ug =Okms', Ver=I, pr=10*gcem3, pr =10-*GPa 


The aluminium is assumed to be purely elastic and is governed by the equation of state from [13] with parameters 
given in Table 1. The air is governed by an ideal gas equation of state with y = 1.4. The problem is run to time 0.6 jus 
with CFL number 0.1. 

Results are shown in Figure 5. The exact solution has been computed for the case with vacuum instead of air, but 
this makes little difference to the solution within the metal which is shown. The air supports a right-going shock not 
present in the vacuum which can be seen in the normal velocity field. 


4.3. Shell collapse 


The shell collapse problem has been studied by many authors, an early contributor being Verney [54], and is a demand- 
ing test of energy conservation and symmetry preservation. This can be seen in solutions obtained using alternative 
Eulerian sharp interface treatments such as the ghost-fluid method, where high resolution is typically required to ob- 
tain a reasonable solution [6, 32]. In contrast our method exhibits smaller conservation errors and is able to provide 
excellent results even on extremely coarse meshes. 

A cylindrical aluminium shell is modelled with initial inner radius of 0.8 cm and initial outer radius 1.0cm; the 
domain size is [0, 1.1]? cm. An initial velocity distribution is applied that causes the shell to collapse incompressibly. 
Plastic distortional effects convert kinetic energy to internal energy, and eventually the shell stops. Howell and Ball 
present formulae for the stopping radii [20]. In this case we have a final inner radius of 0.3cm and a final outer 
radius of approximately 0.671 cm. The initial radial velocity is given by 0.8u9/r where up = 0.437kms~'. The 
equation of state parameters for the aluminium are given in Table 1, and ideal plasticity is used with constant yield 
stress oy = 0.2976GPa. The remainder of the domain is filled with air governed by an ideal gas equation of state 
with y = 1.4 and po = 1.225- 10~* gcm7?. The problem is run to time 30 us with CFL number 0.6. Figure 6 shows 
snapshots of velocity magnitude within the shell, and the location of the interface, at t = 0, 10, 20 and 30 ys. 

The final inner and outer radii are shown as a function of angle in Figure 7a and 7b respectively for a range of 
mesh sizes. In all cases it can be seen that the radial error is limited to less than 2%, even on a very coarse mesh with 
only 64 cells in each coordinate direction. In contrast, both Barton et al. [6] and Lopez Ortega et al. [32] report errors 
in the region of 35% for similarly sized meshes using ghost-fluid methods. At higher (but still moderate) resolutions, 
the deviations from radial symmetry are reduced. 

Conservation histories for total mass and total energy as shown in Figure 7c and 7d respectively for the same range 
of mesh sizes. Although the volume-of-fluid method is not fully conservative, it can be seen that the error in total 
mass for this problem is less than a fifth of a percent, and the error for total energy is limited to less than five percent. 


4.4. Cylinder impact 


Cylinder impacts provide a strenuous test of elasto-plastic flow and are included here to validate the rate-sensitive 
plasticity treatment. A cylinder strikes a rigid anvil (here modelled using a perfectly reflecting boundary) and deforms 
plastically, flowing radially outwards until all kinetic energy has been converted to internal energy and the cylinder 
stops, with a reduced length. Wilkins and Guinan provide experimental data in [58], and we also compare against 
results given by Barton [5] using the underlying diffuse-interface method. 

The initial length of the cylinder is Lo = 2.347cm and the radius is 0.381 cm. We predict the final length Ly 
for a range of impact velocities in 2D axisymmetric geometry, using a domain of 1.25 x 2.5cm and a mesh size of 
512 x 1024. The cylinder is aluminium 6061-T6 with parameters given in Table 1. The Johnson-Cook rate-sensitive 
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Figure 5: Solid-fiuid Riemann problem. 
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Figure 6: The collapsing shell problem: magnitude of velocity within the shell measured in kms~! (velocity of surrounding air not shown) and 
¢ = 0.5 contour at intervals of 10 us from t = Ops to t = 30s. 
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Figure 7: Final radii and energy conservation for the aluminium shell collapse problem at a range of mesh sizes. (a) shows the predicted inner 
radius as a function of the angle @ € [0, 3], (b) shows the predicted outer radius, (c) shows percentage conservation of total mass and (d) shows 
percentage conservation of total energy. The dashed lines show the analytic radii computed using the formulae from Howell and Ball [20]. The 
volume-of-fluid method produces excellent results even at the coarsest mesh resolution, in stark contrast to previously published ghost-fluid method 
results. Radial symmetry is well preserved and improves as the mesh is refined, and conservation errors are low. 
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Material C1 Co C3 n Xo m  To(K)  Twmet (K) 


Aluminium 6061-T6 0.324 0.114 0.002 0.45 1075 - - = 
Copper 0.09 0.292 0.025 0.31 107° 1.09 288.15 1356 


Table 2: Parameters for the Johnson-Cook rate-sensitive plasticity model. 
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Figure 8: The cylinder impact problem with an impact velocity of 373 ms—!: equivalent plastic strain is shown along with the ¢ = 0.5 contour at 
intervals of 10 us from t = Os to t = 30 ys. 


plasticity model is used, including work-hardening effects, with parameters given in Table 2. Thermal softening is 
disabled here by setting Tne) — 00. The cylinder is immersed in quiescent air governed by an ideal gas equation 
of state with y = 1.4 and po = 1.225- 10-2 gcm~?. All materials are initialised to a pressure of | atm and a CFL 
number of 0.6 is used. The time evolution of the cylinder/air interface and the equivalent plastic strain is shown for 
the 373 ms~! case in Figure 8. 

The final ratios L;/Lo for each impact velocity are compared in Figure 9. The volume-of-fluid method matches 
experiment very well, correctly predicting the non-linear trend. 


4.5. Buried explosive 


The buried explosive problem previously studied in [5] is useful to contrast the behaviour of the volume-of-fluid 
method against the diffuse-interface method for problems involving large deformations and interface breakup. A 
charge buried in soft clay detonates resulting in ejecta and crater formation. 

For 2D axisymmetric simulations the domain is [0,2] x [—2,4]m. For full 3D simulations, one quadrant is 
modelled, and the domain is [0,2] x [—2,2] x [0,2]m. The initial location of the clay/air interface is y = 0. The 
material parameters for the clay are given in Table 1. Rather than using a fixed Griineisen parameter, a power-law 
form is used 


q 
r= (4) (104) 


Po 
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Figure 9: Results for the cylinder impact test case. The ratio of the final length of the cylinder to its initial length, both from experiment and 
predicted using the DI+THINC and VOF methods. 


with q = 1. The air is modelled using an ideal gas equation of state with y = 1.4 and py = 1.225-10-* gcem7?. 


The clay and air are given an initial pressure of 1 atm. The explosive products are initialised in the cylinder centred at 
(0, —13.081, 0) cm with radius 8.763 cm and height 5.842 cm. The phenomenological balloon model is used to model 
the explosive products, sharing the equation of state with the air [28]. The initial product density is 1.63 gcm~*, and 
the initial specific internal energy is 6.0568 kJ g~'. Gravitational and atmospheric effects are not important over the 
short timescales considered and so are not included. 

Results for the 2D axisymmetric case are shown in Figure 10. Snapshots of pressure and a numerical Schlieren 
of density are shown at t = 0.3, 1, 1.8 and 3 ms for the volume-of-fluid method and the underlying diffuse-interface 
method with THINC interface sharpening. The function used to define the Schlieren plots is from [43] 


where we take xk = 8000. At the earlier times, excellent agreement is seen between the two methods in terms of cavity 
dimensions and ground shock propagation. At the later times, the crater dimensions agree well and the volume-of- 
fluid method shows its strengths by resolving fine-scale turbulent structures resulting from the break up of the lofted 
clay, which are suppressed in the diffuse-interface solution. A zoomed in density Schlieren at the final time is shown in 
Figure 11a that highlights these details. Closer examination of the log-scale volume-fraction field in Figure 1 1b shows 
that the interface does not fully break when using the diffuse-interface method, instead unphysical thin ligaments with 
low but non-zero volume fraction develop. 

Results from the 3D simulations are shown in Figure 12 including the contours of the clay/air interface and the 
pressure in the ground. Our findings are analogous to the 2D axisymmetric case, with good agreement between the 
two methods on cavity and crater dimensions, as well as ground shock propagation. The 3D volume-of-fluid method 
exhibits the same advantages as in 2D, showing excellent fidelity in the break-up of the lofted clay. 


5. Application: modelling the BRL 105mm shaped charge 


In this section we apply the method to modelling the BRL 105mm unconfined shaped charge as detailed in [2]. 
Shaped charges use convergent detonation waves in a high-explosive to produce a hyper-velocity metal jet with many 
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Figure 10: Comparison between DI+THINC and VOF for the 2D axisymmetric buried explosive problem. The mesh size is 800 x 2400. Panels 
(a), (b), (c) and (d) show t = 0.3, 1.0, 1.8 and 3 ms respectively, and similarly for panels (e)—(h). In each panel the left half of the image shows 
DI+THINC and the right half shows the corresponding VOF solution. The first row shows a numerical Schlieren of density and the second shows 
pressure measured in GPa. The yellow lines denote the ¢ = 0.5 contour. While the two methods agree well on ground shock and crater dimensions, 
it is clear that the DI+THINC method has suppressed the growth of turbulent instabilities in the wake of lofted ejecta, which are much more 
pronounced in the VOF solution. 
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Figure 11: Close-up of tf = 3 ms for the 2D axisymmetric buried explosion problem. (a) shows the numerical Schlieren of density, showing how 
the VOF method captures fine-scale turbulent flow in the wake of the lofted ejecta. (b) compares the log-scale clay ejecta volume fraction fields 
between DI+THINC and VOF. Unphysical ligament formation can be clearly observed in the DI+THINC solution, the interface never truly breaks. 
This is absent in the VOF solution. 


Comp-B po (g cem~?) A(GPa) B(GPa)  R, Ry Tp Cy (kJ K7!¢7!) 
Reactants 1.717 7.781-10*  —5.031 11.3 1.13 0.8938 2.487 - 10-3 
Products 1.717 5.242 102 7.678 4.2 1.1 0.34 10-3 


Table 3: JWL equation of state parameters for Comp-B, taken from [39]. 


industrial and defence applications. They are also stringent tests of elasto-plastic hydrodynamics methods as they 
typically include multiple materials arranged in a complex configuration and subject to extreme deformation. This 
particular shaped charge has been modelled by other authors [57, 3], allowing us to draw useful comparisons. Note 
that it is not our intention to fine tune the model to match experiment, only to show that the volume-of-fluid method is 
capable of treating this extremely strenuous scenario and producing reasonable results. It is noted that the underlying 
diffuse-interface method with THINC interface sharpening does not work well for this problem as the thin jet becomes 
smeared over a large area. 

The 2D axisymmetric domain is [0,6] x [—6, 40] cm, with the origin chosen to be coincident with the inner liner 
apex. The geometry of the shaped charge is given in [2], but unfortunately not all pertinent details are specified. For 
simplicity we omit the tetryl booster pellet and detonator from our model and ignore the rounded apex of the liner. 
Our modified geometry is shown in Figure 13. 

The liner is copper with parameters given in Table 1. The Johnson-Cook rate-sensitive plasticity model is used, 
including both work-hardening and thermal-softening effects, with parameters given in Table 2. The high-explosive 
(Comp-B) is treated using the reactive burn model from [55], summarised in Section 2.4. The reactants and products 
both obey the JWL equation of state, with parameters given in Table 3. The explosive is assumed to have no strength 
hence G = 0. The parameters for the ignition and growth model are a = 0.01, b = 0.222, c = 0.222, d = 0.666, 
e = 0.0,g = 00, x = 40,y = 2.0,z = 00, Fy = 0.3, Fe, = 1.0, Fe, = 10,7 = 4.4- 107 (10us)—!, 
G, = 4.14- 107! GPa-*(10us)~! and G, = 0.0 GPa ~*(10us)~!, taken from [39]. The heat of detonation associated 
with the Comp-B is 6.3 kJ g~'. The shaped charge is immersed in air governed by an ideal gas equation of state with 
y = 1.4 and pp = 1.225- 10-7 gem7?. 

All materials are initialised to a pressure of | atm, except for a hemispherical booster region centred at the back of 
the explosive with radius 1.439 cm and pressure 27 GPa, used to initiate the detonation. Twenty five Lagrangian tracer 
particles are initially located along the inner surface of the liner, evenly spaced along the z-axis at intervals of 0.25cm 
starting from z = 3.0cm. The time history of the velocity is measured for each tracer particle. 

Figure 14a—14d shows the high-explosive burn and subsequent collapse of the liner that occurs during the first 
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Figure 12: Comparison between DI+THINC and VOF for the 3D buried explosive problem. The mesh size is 240 x 480 x 240. (a)-(d) show 
t = 0.3, 1.0, 1.4 and 1.8 ms respectively. In each image the left half shows DI+THINC and the right half shows the corresponding VOF solution. 
The yellow contour shows the clay/air interface and the colouring indicates pressure. As in 2D the two methods agree well on ground shock and 
crater dimensions, but VOF offers much improved fidelity for the lofted ejecta. 


© British Crown Owned Copyright 2023/AWE 29 


6.0 in 


3.4in x 3.25 in 


Comp-B Copper 


Figure 13: The BRL 105mm unconfined shaped charge, simplified from the diagram given in [2]. 


eighteen microseconds. The detonation wave is initiated by the high pressure booster region. Just after f = 3 ys the 
wave strikes the apex of the liner and pressure begins to build up in that region. At ¢ = 12s the jet is beginning 
to form as the liner turns in and is rapidly accelerated. By t = 18 ys the high-explosive is fully burnt and the high 
pressure stagnant region can be clearly seen. Figure 14e and 14f show the late-time jet formation at t = 45 and 
60 ys respectively. The jet reaches a stable speed and begins to thin out due to the velocity gradient along its length. 
The jet also has a blunt tip, which is in good agreement with the Arbitrary Lagrangian-Eulerian results presented by 
Barlow et al. [3], and as they note matches what is typically seen in experiment. Close inspection of the density field 
reveals fragments of copper breaking off the tip. 

Figure 15a shows the time history of the velocity of each Lagrangian tracer particle over the course of the simu- 
lation. The inner apex of the liner is accelerated close to the detonation velocity of the Comp-B (7.98 mm pus” !). By 
the end of the simulation the particle velocities have stabilised. Figure 15b shows the slice corresponding to t = 60 ys 
as a function of the initial z-coordinate of each tracer particle. The experimental results from [2] are overlaid. The 
jet velocity is particularly sensitive to the parameters used in the reactive burn model which have not been tuned, but 
nevertheless a fairly good match is obtained. 


6. Conclusions 


We have presented a novel Godunoy-type volume-of-fluid method for the simulation of large deformation multi- 
material problems featuring elasto-plastic solids and fluids. The method can be viewed as an interface-sharpening 
extension to the diffuse-interface method of Barton [5]. The method is practical, robust and extensible, and has been 
shown capable of treating a range of challenging problems with excellent fidelity. 


References 


[1] G. Allaire, S. Clerc, and S. Kokh. A five-equation model for the simulation of interfaces between compressible fluids. Journal of Computa- 

tional Physics, 181(2):577-616, 2002. 

[2] FE. Allison and R. Vitali. An application of the jet-formation theory to a 105mm shaped charge. Technical Report AD0277458, Army 

Ballistic Research Lab, Aberdeen Proving Ground, MD, 1962. 

[3] A.J. Barlow, R.N. Hill, and M.J. Shashkov. Constrained optimization framework for interface-aware sub-scale dynamics closure model for 

multimaterial cells in lagrangian and arbitrary lagrangian—eulerian hydrodynamics. Journal of Computational Physics, 276:92-135, 2014. 

[4] A.J. Barlow, P.H. Maire, WJ. Rider, R.N. Rieben, and M.J. Shashkov. Arbitrary lagrangian—eulerian methods for modeling high-speed 

compressible multimaterial flows. Journal of Computational Physics, 322:603-665, 2016. 

[5] P.T. Barton. An interface-capturing godunov method for the simulation of compressible solid-fluid problems. Journal of Computational 
Physics, 390:25—50, 2019. 


© British Crown Owned Copyright 2023/AWE 30 


(f) 
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In each panel, the top half shows pressure measured in GPa and the bottom half shows density measured in g cm~*. (e) and (f) show the late-time 
jet formation at t = 45 and 60s. The top half of each panel shows the axial speed measured in mmps~! and the bottom half shows the density. 
The yellow line denotes the @ = 0.5 contour. 
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